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I.  INTRODUCTION 


Factorial  series  for  Bessel  functions,  confluent  hypergeometric 

functions,  and  certain  other  special  functions  can  be  used  to  check  the 

accuracy  of  calculations  for  these  functions  provided  the  argument  is  not  too 

small.  Generally,  if  a  function  is  analytic  in  the  right  half  plane, 

including  the  imaginary  axis,  and  can  be  represented  as  a  Laplace  transform, 

then  a  factorial  series  can  be  derived  which  will  converge  in  the  right  half 

plane.  Buchal  and  Duffy*  obtained  factorial  series  for  Hankel  functions. 

Coulomb  wave  functions,  and  Mathieu  functions.  These  authors  also  studied  the 

convergence  properties  of  the  factorial  series  for  the  Hankel  functions  in 

considerable  detail.  Their  results  showed  that  factorial  series  were  more 

accurate  than  Hankel  asymptotic  series  for  the  same  argument.  The  analysis 

o 

was  based  on  Bernoulli  polynomials  as  discussed  by  Doetsch  and  Milne- 
3 

Thomson  . 

The  analysis  in  this  paper  is  based  on  an  algorithm  of  Wasow^  which  uses 
Stirling  numbers  of  the  first  kind  to  calculate  coefficients  for  the  factorial 
series.  The  Stirling  numbers  of  the  first  kind  increase  very  rapidly  with 
order,  eventually  obtaining  overrun  in  the  electronic  computer.  Moreover,  the 
factorial  series  for  complex  argument  is  awkward  if  there  is  a  branch  point  at 
the  origin.  Rosser'’  obtained  a  generalized  factorial  series  for  modified 


R.N.  Buchal  and.  G.  Duffy,  "Factorial  Series  Representations  as  Check 
Solutions  in  the  Development  of  Special  Function  Subroutines,"  Argonne 
Rational  Laboratory  Technical  Memorandum  No.  200,  May  1970. 

2 

G.  Doetsch,  Hanbuch  Per  Laplace-Trans formation.  Band  II,  Birkhauser  Verlag, 
Rosel  and  Stuggart,  1955,  pages  201-232. 


3 

L.M.  Milne  -  Thomson,  The  Calculus  of  Finite  Differences.  Macmillan  and  Co., 
Limited,  London,  1933. 


4W.  Wasow,  Asymptotic  Expansions  for  Ordinary  Differential  Equations . 
Interscience  Publishers,  a  Division  of  John  Wiley  &  Sons,  Inc.,  New  York, 
pages  323-331. 


5-/.  3.  Rosser,  "Factorial  Expansions  for  Certain  Bessel  Functions," 
Mathematics  Research  Center  Summary  Report  No.  1572,  November  1975. 
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Bessel  functions  of  the  second  kind  by  direct  manipulation  of  the  Laplace 
transform  and  also  established  the  convergence  properties.  His  analysis  is 
quite  difficult  and  requires  a  separate  treatment  for  each  case.  In  this 
paper  we  derive  Stirling  numbers  of  the  first  kind  and  fractional  order, 
leading  to  a  generalization  of  Wasow's  algorithm.  Overrun  in  the  computer  is 
eliminated  by  scaling  the  Stirling  numbers  of  integral  and  fractional  order  by 
omitting  the  gamma  functions  in  the  definition. 

A  further  generalization  of  factorial  series  is  required  if  the  Laplace 
integral  representing  the  function  is  evaluated  between  finite  limits.  If  a 
factorial  series  is  regarded  as  a  series  of  beta  functions,  it  is  logical  to 
represent  a  Laplace  integral  with  finite  limits  in  terms  of  a  series  of 
incomplete  beta  functions.  The  new  series  obtained  by  this  method  is 
convergent  even  though  the  corresponding  factorial  series  may  be  divergent. 

II.  MODIFIED  STIRLING  NUMBERS  OF  THE  FIRST  KIND  AND  FRACTIONAL  ORDER 

Stirling  numbers  of  the  first  kind  are  defined  as  coefficients  which 
occur  when  a  factorial  is  expanded  into  a  polynomial^ *^ ,  given  by 

x  (x-1)  (x-2)  ...  (x-n+1)  ■  \  S^™1^  h™  •  (1) 

m=o 


Clearly  m  must  be  an  integer  in  the  above  equation.  However,  a  generating 
function  involving  logarithms  is  not  subject  to  this  restriction,  such  as 


00  .  *  n 

[in  (1+x) ]m  *  m!  \  Sn  m  ~  ,  for  1xf<l 
n-m 


(2) 


fi  # 

M.  Abramowitz  and  L.A.  Stegun,  Handbook  of  Mathematical  Functions,  Graphs, 
and  Mathematical  Tables .  National  Bureau  of  Standards,  Washington,  D.C., 
No.  56,  Applied  Mathematics  Series,  June  1964.  See  Formulas  24.1.3,  page 
824,  for  discussion  of  Stirling  numbers .  Continued  fractions  for  the 
incomplete  beta  functions  are  given  in  Formulas  26.5.8,  page  944. 

n 

C.  Jordan,  Calculus  of  Finite  Differences,  Rottig  and  Romwalter ,  Sapron, 
Hungary,  1939.  Stirling  numbers  are  discussed  in  Chapter  4,  pages  142-229. 


or,  alternatively, 


Ita  ci«)l»  -  r(«i)  £  sn<“> 

n=m 

In  order  to  avoid  the  gamma  function  and  to  include  non-integral  indices  we 
define  W  ^  ^  by  the  equation 


[in  (l+x)]U  -  l  Wk(u)  x 
k-o 


(3) 


where  u  may  take  on  fractional  as  well  as  integral  values.  When  the 

have  been  calculated,  Stirling  numbers  of  the  first  kind  may  be  obtained  from 

the  following  equations. 


S  (m)  -  r  (n+1)  W<“>  /  r(m+l)  , 

n  n-m 


(4) 


where 


u  =  m  (5) 

k  =  n  -  m  .  (6) 


To  find  *  divide  Eq.  (3)  by  xV  and  evaluate  the  limit  of  each  side 

of  resulting  equations  as  x  ♦  0.  We  obtain 


(7) 
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We  can  also  prove  that 


0,  for  k  >  1 


(8) 


and 


(9) 


We  now  derive  a 

,(u)  w(°) 

•  •  •  *  » 


sequence  of  triangular  equations  for  calculating 
In  turn.  Let 


y(x)  *■  [An  (l+x)J 


(10) 


and 


u(x)  «  (x+1)  An  (1+x). 


(11) 


Then 


u(x)  y'(x)  »  oy(x)  . 

i 


(12) 


But 


y(x) 


l 

k“o 


u+k 

x 


(13) 


y'(x) 


l  „<u>  fu+k)  X  u+k_1  and 
L  k 


(14) 
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for  k>l 


(15) 


\  ,  1  /  2  13,1  4 

u(x)  -  X+V2X  -  -r  x  +y2x  (“1 ) 


k(k-l)  » 


On  inserting  these  series  into  Eq.  (12),  carrying  out  the  indicated 
multiplication  and  equating  coefficients  of  like  powers  of  x  on  each  side  of 
the  resulting  equation,  we  find 


w(u)  *  w(u), 

o  o  ' 


^  w<u>  +  l/2  O  -  U  W<U)  and 


w(y)  +  «±L  w(y)  .  «  .  „  w(«) 

1  2  2  wl  6  w2 


On  rearranging  these  equations  we  find 


w[U  +  V2  U  -  o  , 


-  o. 


3  w<")  +  »«  „<">  -  „<”>  + 

J  l  i  0  1 


and,  in  general, 


,  \  *-l 

*W(U>  +  I  (-1) 
k-0 


i+k+1  v+k 


(A-k)(i-k+l)  k 
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r 

f 


To  find  a  recurrence  formula  involving  different  orders,  differentiate 
Eq.  (3)  with  respect  to  x  and  then  multiply  both  sides  of  the  resulting 
equation  by  (x+1),  giving 


u  [Jin  (l+x)]U_1  =  l  w£U)  (u+k)(x+l)  x  U+k  1  .  (20) 

k=o 


If  we  replace  u  by  (u-1)  in  Eq.  (3)  and  multiply  both  sides  of  the  resulting 
equation  by  u  ,  we  find 


u  ( An( 1+x)] 


u-1 


l 

k=o 


„<^D  0][Ufk-i 
k 


(21) 


On  comparing  coefficients  of  like  powers  of  x  in  the  last  two  equations  we 
find 


/U) 

k+i 


=  [o  w, 


(u-1) 

k+1 


-  (u+k)  W<U)]  /  [u+x+1]. 


(22) 


Hence  we  can  calculate 


fro. 

k+1 


W(U)  and  W(UH) 
k+1  and  Wk 


Higher  order  numbers 


can  be  generated  in  succession  in  the  same  manner.  The  values  of 

must  be  calculated  from  Eq.  (19)  before  the  recurrence 
formula  given  by  Eq.  (22)  can  be  used;  a  double  entry  table  is  finally 
obtained . 


Wasow  uses  Schloralich's  definition  of  factorial  coefficients  in  his 
development  of  factorial  series  ^ 


n-1 

x  (x+1)  (x+2)  ...  (x+n-1)  =  l  T  n  xn  m  . 

m 

ra=o 


(23) 


It  follows  that  the  factorial  coefficients  and  Stirling  numbers  of  the  first 
kind  are  related  by  the  formula 


r  n  -  (-i)  n_m  s<m) 

n-m  n 


We  define 


V(m)  =  (-l)m"n  r(nri-l)  s(m)  /  r(n+i ) 

n-m  n 


and  for  fractional  orders 


,  /,  x.u  „(u)  u+k 

in(l-x)]  ■  i  V.  x  . 


Finally,  we  obtain  the  following  recurrence  formulas  in  the  manner  indicated 
previously 


i  V(u)  -  V  — ** _ V(u)  -  0  and 

1  k-o  k 


Vk+1  "  lu  Vk+1°  +  (o+k)  VkU)  1  1  tu+k+1)  • 


III.  GENERALIZED  FACTORIAL  SERIES 


We  now  derive  an  extension  of  Wasow's  algorithm  for  a  Laplace  integral  to 
functions  with  a  branch  point  at  the  origin.  The  branch  point  involves 

Q 

fractional  powers,  in  the  same  context  as  Watson's  Lemma  ;  logarithmic  branch 
points  are  not  considered.  Assume 


°E.T.  Copeon,  An  Introduetion  to  the  2j 
Variable .  Oxford,  1935,  pages  218-219. 
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F(x)  -  /  f(t)  e_Xt  dt  (29) 

o 


and  let 


t  -  -  in  (1-u)  ; 


(30) 


then 


F(x)  =  /  f[-in(l-u)]  [l-u]X  *du 


(31) 


If 


f(t)  »  y  at  U+n,  for  u>-l  and  0<t<l  , 
n 

n**o 


(32) 


then 


<30  \ 

F(x)  =*  I  a  /  [-in(l-u)]  U+n[l-u]X  ^u. 

n 

n=o  o 


(33) 


On  referring  to  Eq.  (26)  we  see 


00  00 


F(x) 


r  v  „(u+n)  ,  u+k+n  s  x-1, 

I  i  an  Vk  ■*  u  (1_u)  du* 


n-o  k-o 


(34) 


Now 


14 


B  (a, 8) 


(35) 


so  that 


t  0-1  (l-t)B_1dt 


F(x)  =  l  l  an  VkU+n)  B  (k+n+U+l,x) 


n=o  k=*o 


(36) 


We  enter  the  terms  in  a  double  entry  table;  then  by  summing  the  diagonals,  we 

q 

obtain  the  Cauchy  sum  of  the  double  series  .  Let 


k+n=£ 


(37) 


and  define 


b 


a 


i 

l 


n-o 


a 

n 


v(u+n)  # 
*-n 


(38) 


then 


F(x)  ■»  l  b  B(u+A+l,x)  . 
S.-0 


(39) 


On  noting  that 


B(o+£+l,x) 


B  (u.x) 


u  (u+1)  ...  (u+&) _ 

(x+u)(x+u+l) . .  .(x+u+Jl) 


(40) 


9K.  Knapp,  Infinite  Sequences  and  Series,  Dover  Publications,  New  York,  1956, 
pages  82-85. 
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and  using  Pochhammer's  symbol  to  represent  the  factorials,  we  find 


F(x)  -  B  (u,x)  l  b  <v>)  .  /  (x+u). 


in  a  formal  sense.  By  analogy  with  conventional  factorial  series,  the  series 
should  converge  in  a  half  plane  which  lies  to  the  right  of  the  imaginary 
axis.  Details  of  the  required  analysis  will  not  be  considered  at  this  time. 

IV.  A  SERIES  OF  INCOMPLETE  BETA  FUNCTIONS 

Since  a  generalized  factorial  series  is  in  fact  a  series  of  beta 
functions,  it  is  natural  to  represent  a  Laplace  integral  with  finite  limits  of 
integration  as  a  series  of  incomplete  beta  functions.  The  lower  limit  of 
integration  can  be  taken  equal  to  zero  without  loss  of  generality.  Assume 


F(x)  *  /  e  tX  f(t)  dt,  t>0  . 


.  -t 

e  ■  1-e  ; 


F(x)  -  /  fl-An(l-u)]  [l-u]X-1du 


On  referring  to  Eq.  (26)  we  find 
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F(,>  -  2  l 

n“o  k«o 


r  u+k+n,,  vX-1 , 

J  u  ( 1-u)  du 


(45) 


Since 


B  (a,B)  =  /  t  a  1  (1-t)  ^  1  dt  , 


(46) 


we  find 


00  00 


F(x)  =  \  1  a  V,  ^V+n^B  (u+k+n+l,x) 


n=o  k=o 


n  k  e 


(47) 


or 


f(x)  =  [  bt  Be(v+«,+l,x) 


l=o 


(48) 


on  referring  to  Eqs.  (37)  and  (38). 

To  compute  the  incomplete  beta  function,  set 


\H-Jl  -  6 


(49) 


and  use  the  formula 


B  (6+1, x)  -  B(x  ,6+1 )  -  B,  (x ,6+1 )  . 
£  1 


(50) 


The  beta  function  on  the  right  side  of  Eq.  (50)  was  expressed  in  terms  of 
gamma  functions,  as  shown  in  Eq.  (51), 
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b(x,6+d  -  r(x)  r(6+i)  /  r(x+6+n 


(51) 


The  gamma  functions  were  obtained  from  the  subroutine  CDLGAM  by  H.  Kuki*^. 

This  subroutine  Is  valid  for  both  real  and  complex  values  of  x. 

The  incomplete  beta  function  on  the  right  side  of  Eq.  (50)  is  given  by 

the  Integral  formula 


B^U,  5+1) 


1-e  . 

r  X- 1  ,  ,  .6, 

J  t  ( 1 — t )  dt. 


(52) 


On  expanding  the  binomial  factor  in  the  Integrand  and  Integrating  term  by  term 
we  find 


B1_e(x«6+1)  55 


(1-e)7 


6_ 

1! 


(1-e) 

(x+1) 


x+1 


6(6-1)  (1-e) 


x+2 


2! 


(x+2) 


6(6-l)(6-2) 

3! 


(1-e) 

(x+3) 


x+3 


+., 


(53) 


This  series  is  satisfactory  if  6  Is  small,  but  is  subject  to  round-off  error 
if  6  is  large  and  positive.  To  overcome  this  difficulty,  integrate  the  right- 
hand  side  of  Eq.  (52)  repeatedly  by  parts.  We  find 


B1_e(x,6+1) 


( 1-e ) 


6  6(l-e) 

x ( x+ 1 ) 


x+1 


6-1 


+  * (6-1) (1-e)**2 

x(x+l)(x+2) 


6-2  .  .  6(6-l)(6-2)...(6-m)  _ 

’**  x(x+l )(x+2) . . . (x+m)  a 


(54) 


where 


10H.  Kuki,  Communiaationa  of  the  ACM.  Vol  IS,  No  4,  April  1972.  The 
subroutine  CDLGAM  i Jae  extracted  from  Algorithm  421. 
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(55) 


R  -  B,  (x+m+1 ,  6-m) 
m  I- e 


We  choose  m  so  that 


1  <  6-m  <  2.  (5b) 

Then  the  terms  of  the  series  given  by  Eq.  (54)  are  positive  when  x  is  real  and 
positive,  and  consequently  the  round-off  error  should  be  small. 

The  remainder  R  is  calculated  from  Eq.  (53), 


( l-e)x+m+*  (  <S-m-l )  ( 1-e  )x+m't'^  (6-m-l)(6-m-2)(  1-e)  +m 

m  x+m+1  1!  (x+m+2)  2!  (x+m+3) 


All  the  terms  of  this  series  after  the  first  term  are  negative,  and  decrease 
rapidly  in  magnitude.  Hence,  the  error  in  calculating  Rm  should  be  small. 

V.  CALCULATIONS 

We  calculated  the  modified  Bessel  functions  KQ(x),  K^(x),  IQ(x),  and 
I ^ (x )  from  the  integrals 


r(  V2 )  e  x  ~ 
K  (x)  - - - -  I 

°  r(  V2  )  o 


/  e~Xt(t2+2t)' L/2dt  , 


KL(x) 


r(V2)  (l/2*)« 


/  e-xt(t2+2t)  ‘/2dt. 


I  (x)  .  - j2  e"xt(2t-t2)'i/2dt, 

°  r(V2)r(V2)  o 


ir(x) 


(  V2x)  eX  2  ,  I, 

- _  /  e  Xt(2t-t2)  '2  dt . 

r(  V2 )  r<-)  O 

£  2 


IS 


On  referring  to  Eqs.  (29)  and  (32),  we  see  that 


(62) 


in  Eqs.  (58)  and  (60),  and 

v  =  V2  (63) 

in  Eqs.  (59)  and  (61).  Hence,  the  modified  Stirling  numbers  V  ^2 +ra  an<j 
Vk2  m  are  required  for  the  coefficients  of  the  factorial  series.  A  short 
table  of  these  numbers  is  given  in  Table  1 . 

Next,  the  coefficients  an  for  the  series  expansions  of  (1+ V2  O  ^2  » 

(  1-  V2  t)~  ^2  t  (  1+  V2  t)  ^2  and  ( 1-  V2  t)  ^2  were  calculated  from  the  appropriate 
recurrence  formulas,  as  shown  in  Table  2. 

The  coefficients  b^  were  calculated  from  Eq.  (38)  for  each  of  the  four 
cases  listed  above.  In  addition,  the  partial  sums 


C 

m 


m 


(64) 


were  calculated  in  order  to  study  the  convergence  of  the  factorial  series. 

00 

The  series  £  Cm  must  converge  if  the  corresponding  factorial  series  is  to 
m 

converge.  This  condition  is  apparently  violated  for  the  functions  IQ(x)  and 
Ij(x),  which  shows  why  the  factorial  series  for  these  functions  apparently 
diverged.  These  results  are  given  in  Table  3. 

The  beta  functions  and  incomplete  beta  functions  required  in  the  series 
expansions  for  KQ(x),  Kj(x),  IQ(x),  and  Ij(x)  were  calculated  from  formulas 
discussed  previously.  Sample  tabulations  are  shown  in  Table  4.  Finally,  the 
modified  Bessel  functions  were  calculated  for  a  limited  range  of  variables. 
These  results  are  shown  in  Table  5. 
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VI.  RESULTS  AND  CONCLUSIONS 


These  series  expansions  in  terms  of  beta  functions  and  incomplete  beta 
functions  were  derived  in  order  to  check  the  accuracy  of  our  Bessel  function 
subroutine  ’  with  independent  calculations.  The  error  analysis  of  our 
subroutine  by  theoretical  methods  would  be  very  difficult,  especially  for  the 
section  involving  continued  fractions.  Hence,  computational  efficiency  is  a 
secondary  consideration  for  the  new  series  expansions.  Addressing  the  comment 
of  the  reviewer*,  we  believe  our  algorithm  for  the  incomplete  beta  function  is 
as  efficient  as  the  continued  fractions  of  Stegun^  when  x  is  large,  since  Eq. 
(54)  can  then  be  used  without  the  remainder.  We  have  not  made  any  specific 
comparison  for  small  values  of  x. 

Since  factorial  series  are  a  method  of  summing  certain  asymptotic  series, 
they  are  most  effective  for  large  and  moderately  large  values  of  the 
argument.  The  convergence  is  slow  when  x  is  small,  so  that  an  excessive 
number  of  terras  is  required.  Round-off  error  may  occur  in  the  coefficients  b^ 
and  in  summing  the  series.  Alternate  methods  of  calculating  Bessel  functions, 
such  as  quadratures,  are  required  when  x  is  small.  Subroutines  used  for 
checking  should  not  use  continued  fractions  or  other  procedures  used  in  the 
subroutine,  to  insure  the  calculations  are  in  fact  independent. 

The  generalized  factorial  series  for  KQ(x)  and  Kj(x)  yield  accurate 
numerical  values  when  x  is  only  moderately  large  and  the  Hankel  asymptotic 
series  is  not  sufficiently  accurate.  On  the  other  hand,  the  new  series  for 
I0(x)  and  Ij(x)  have  the  same  range  of  accuracy  as  the  Hankel  asymptotic 
series,  and  do  not  offer  any  computational  advantage.  This  is  probably  due  to 
the  apparent  divergence  of  the  series  for  the  partial  sums  of  b^. 

The  generalized  factorial  series  for  KQ(x)  and  K^(x)  are  being  extended 
to  the  complex  plane.  Programming  of  generalized  factorial  series  for  the 
ordinary  Bessel  functions  is  in  progress.  Alternate  methods  of  calculating 


H A.  S.  Elder,  "Formulas  for  Calculating  Bessel  Functions  of  Integral  Order 
and  Complex  Argument, "  BEL  Report  No.  1423,  November,  1968. 
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In(x)  are  also  being  considered,  as  the  results  obtained  in  this  paper  fell 
short  of  our  expectations. 

The  modified  Stirling  numbers,  as  defined  in  this  paper,  are  more  useful 

for  computations  involving  Wasow's  algorithm  than  the  original  Stirling 

numbers,  as  problems  arising  from  very  large  numbers  and  overrun  of  the 

computer  registers  are  entirely  avoided.  The  method  of  scaling  employed  in 

this  paper,  which  merely  involves  the  omission  of  gamma  functions  in  the 

definition  of  Stirling  numbers,  is  more  effective  than  the  method  of  scaling 

1 3 

used  in  previous  paper  by  the  authors 
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